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I. INTRODUCTION 



Strongly interacting matter at non-zero baryon density and high temperature is a subject of great interest for a 
wide spectrum of physicists. A deep understanding of the different facets of strongly interacting matter, especially, 
the physics of colour deconfinement might help us to get a better picture of various astrophysical and cosmological 
phenomena. Huge accelerators have been built at CERN, Geneva Q and at RHIC, Brookhaven |2| to collide beams 
of heavy nuclei at relativistic energies to recreate the matter in such extreme conditions. The focus of the heavy- ion 
physics program is to study and understand properties like the emergence of macroscopic collective phenomenon 
and their role on the evolution of a strongly interacting system. The analysis of data so obtained requires a proper 
understanding of Quantum Chromodynamics (QCD) , the theory of strong interactions. The thermodynamic properties 
are obtained from the QCD equation of state (EOS) and various transport coefficients. 

The weak coupling expansion of the free energy within perturbative QCD (pQCD), is presently known Q to order 

as or g . However, in spite of the higher order, the series show a deceptively poor convergence except for coupling 
constants as low as as < 0.05, corresponding to temperature as high as 10^ Tc. 

There has been significant recent activity to improve on this convergence by Hard Thermal Loop (HTL) re- 
summation schemes |^ which give a description of the plasma in terms of weakly interacting hard and soft quasi- 
particles 0, 0, • The former are massive excitations with masses mo ^ gT, while latter are either on-shell collective 
excitations or virtual quanta exchange in the interactions between hard particles. In contrast to ordinary perturbation 
theory, the HTL and its next-to- leading order approximations are well controlled in the regime above 5Tc, which is how- 
ever, much higher than the temperature range in present day heavy-ion collisions. In addition the different approaches 
within the same general framework of HTL approximations lead to different predictions [8Ll9Lll0llllL IT^IT3Lll4| . and 
cannot be made fully systematic when applied to thermodynamics. At the same time, this approach lacks a reliable 
and comprehensive description below 5Tc due to their stronger sensitivity to non-perturbative collective phenomena. 

There is also a non-perturbative first principle method to compute the contribution of the soft fields to the thermo- 
dynamics. At high temperatures the hard modes get large temperature dependent masses, which are then integrated 
out leaving a three-dimensional effective theory. This method is known as dimensional reduction jl5l Il6l ITtI I . How- 
ever, as expected, this method is not well-suited for lower temperature and the susceptibilities calculated from these 
method [ig do not agree with those of Lattice QCD computations below ~ 2Tc. 

Exploring the qualitative features of strongly interacting matter and making quantitative predictions of its properties 
from first principles is the central goal of numerical studies of equilibrium thermodynamics of QCD within the 
framework of Lattice regularization. Over the years this formulation has given us a wealth of information (for a 
review see Ref.^^)- It has now been established that there is only a crossover of normal hadronic matter to a 
state of deconfined quarks and gluons at temperature Tc ~ 200 MeV. At this temperature, deconfinement of color 
charge and restoration of chiral symmetry is found to occur simultaneously. Moreover, the equation of state, various 
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susceptibilities and transport coefficients have been obtained. To perform Lattice QCD computations at non-zero 
cliemical potential however is still a non-trivial task due to the complex fermion determinant. Recently thou gh, there 
have been various new developments to tackle this problem for small chemical potentials |20, UM, 12^ l2M ■ 
Generally the QCD inspired phenomenological models are much easier to handle compared to Lattice or the per- 
turbative QCD calculations. But in all these models despite their simplicity, the absence of a proper order parameter 
for dcconfincmcnt transition adds to the uncertainties inherent in such studies and hence reduces the predictive power 
of such models. 

The thermal average of the Polyakov loop can be considered to be the order parameter for deconfinement transition 
[27j . Hence a judicious use of the Polyakov loop in effective models may prove to be of great advantage. Some 
important developments have taken place in obtaining effective theories for Polyakov loop from that of the temporal 
background gauge field Isof . It is expected that such an effective theory for the Polyakov loop should work 

at least near the phase transition where long distance physics become important. More recently, the parameters in 
these effective theories have been fixed |^ using the Lattice data (similar comparisons of perturbative effects on 
Polyakov loop with Lattice data above the deconfinement transition was studied in |23)- These Polyakov loop models 
have been applied in the context of cosmology . 

For non-zero chemical potential there are various QCD inspired models which indicate (see e.g. Refs. [ssl Is^ IbtI IsM 
Isgf ) that at low temperatures there is a possibility of first order phase transition for a large baryon chemical potential 
fj,Bc- This fiB^ is supposed to decrease with increasing temperature. Thus there is a first order phase transition 
line starting from {T = 0, /is = /is^) on the fj,B axis in the {T,^b) phase diagram which steadily bends towards 
the {T = Tc, hb = 0) point and may actually terminate at a critical end point (CEP) characterized by (T = Tg, 
Ms — I^Be)^ which can be detected via enhanced critical fluctuations in heavy-ion reactions ^^l- The location of this 
CEP has become a topic of major importance in effective model studies (see e.g. Ref.^J). Also on the Lattice the 
CEP was located for the physical [i^ and for somewhat larger quark masses using the reweighting technique of 
[23 |. and for Taylor expansion method in [25l |. 

Thus the broad picture of QCD thermodynamics is becoming more transparent with the combination of research 
in the areas of perturbative QCD, effective models and Lattice computations. In this paper we study some of the 
thermodynamic properties of strongly interacting matter using the Polyakov loop -I- Nambu-Jona-Lasinio (PNJL) 
model |4j. The motivation behind the PNJL model is to couple the chiral and deconfinement order parameters 
inside a single framework. Similar approach to understand QCD thermodynamics is being actively pursued by various 
authors (see e.g. Ref. 45] and references therein). Here we shall first compute the EOS and the quark number 
susceptibilities. Susceptibilities in general, are related to fiuctuations via the fiuctuation-dissipation theorem. The 
difference in fluctuations of various conserved quantities like baryon number, electric charge etc. in the hadronic and 
deconfined phases is supposed to signal the phase transition between these two phases in heavy- ion reactions |46| . 
Measurements of these fluctuations have taken a central place in the heavy-ion collisions 1471. Recently, computations 
on the Lattice have obtained many of these susceptibilities at zero chemical potential |25l |2^ . and our flrst aim 
would be to compare the corresponding quantities extracted from the PNJL model. As we shall show here that the 
agreement of the values obtained from the PNJL model is quite satisfactory with the Lattice data. 

Subsequently we have studied few more quantities of interest. One such quantity is the specific heat Cy, which 
is related to the event-by-event temperature fluctuations 48], and mean transverse momentum fiuctuations |49j in 
heavy-ion reactions. These fluctuations should show a diverging behaviour near the CEP. Next, we obtain the speed 
of sound (basically its square, w^), which determines the flow properties in heavy-ion reactions .50, 51, 52, 5^]. Using 
the proper hydrodynamic equations including the speed of sound it is possible to analyze the rapidity distribution of 
secondary particles in collision experiments |54|. Finally we obtain the conformal measure C = A/e, where A = e — 3P 
is the interaction measure and e and P are respectively the energy density and pressure of strongly interacting matter. 
As has been pointed out in [55l Is^ , the conformal measure seems to be emerging as an important measure to draw 
similarities between long distance physics of QCD and conformal field theory, with results coming from both the areas 
of AdS/CFT correspondence [H^], RHIC data [E^ and Lattice computations fsoj . 

The plan of the paper is as follows. In section 2, we discuss the PNJL model briefly. In section 3, we present 
the formalisms. The EOS has been obtained from the Lagrangian of PNJL model in mcan-flcld approximation. The 
Taylor expansion coefficients of pressure with respect to the quark number chemical potential /io, is written down. 
We also present the various formulae for Cy, Vg and C. In section 4, we present our results and comparisons with 
Lattice data. Finally we conclude with a discussion in section 5. 



II. PNJL MODEL 



The PNJL model was first introduced in Ref.'S] to couple the Nambu-Jona-Lasinio (NJL) mode l IgoI] (see Ref.|61| 
for most recent developments), with the Polyakov loop. Recently, this has been extended in Refs.|63.l63| to include the 
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Polyakov loop eflFective potential |^ [s^l ■ While the NJL part is supposed to give the correct chiral properties, 

the Polyakov loop part should simulate the deconfinement physics. Indeed such a synthesis worked well to firstly 
confirm the "coincidence" of onset of chiral restoration and deconfinement as observed in Lattice simulations (see 
discussions in Kef's?!) in the PNJL model ^6^. Secondly, the pressure, scaled pressure difference, number density 
and the interaction measure were computed in Ref. 63J for two quark flavours, and all the quantities compared well 
with the lattice data. This model is supposed to work upto an upper limit of temperature, because here the gluon 
physics is contained only in a static background field that comes in the Polyakov loop. However, transverse degrees 
of freedom will be important for T > 2.5Tc |65l |. 

For the details of the PNJL model parameterization, the reader is referred to Ref. j63| . Our starting point is the 
thermodynamic potential per unit volume given by. 



n = U {<i>,^,T) + — - 2NfT 



(2^ 



I In l + 3(c 



In 
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(1) 



Here, U ($,$,T) is the effective potential for the traced Polyakov loop <i> and its conjugate <&, and T is the 
temperature. The functional form of the potential is, 



W($,4',r) 
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with 



&2(T) = ao + a,(|)+a.(|) + as (|) ■ (3) 

The coefficients ai and hi were fitted from Lattice data of pure gauge theory. The parameter Tq is precisely the 
transition temperature for this theory, and as indicated by Lattice data its value was chosen to be 270 MeV 'b^.'btIIb^. 
With the coupling to NJL model the transition doesn't remain first order. In this case from the peak in d^/dT the 
transition (or crossover) temperature comes around 230 MeV. The authors in Ref.ls^ then reduce the To to 190 
MeV such that the Tc becomes about 180 MeV, commensurate with Lattice data with two flavours of dynamical 
fermions 69]. However, we shall keep using Tg = 270 MeV, since for Tq = 190 MeV, there is about 25 MeV shift in 
the chiral and deconfinement transitions with all other model parameters remaining fixed, as compared to less than 5 
MeV shift for Tq = 270 MeV. We have checked that for both the values of Tq, the susceptibilities we measure, when 
plotted against T /T^. show very little dependence on Tq. 

The other notations in are as follows, a is the auxiliary field introduced via bosonization techniques in the 
PNJL Lagrangian. < a >= G < ipip > is the chiral condensate. G is the effective coupling strength of a local, 
chiral symmetric four-point interaction. Nf denotes number of flavours. Ep ~ yTT^ + m^, where m = mo — (cr) = 
mo — G{ipTp), with mo = m„ = rrid being the value of current quark mass. A is the 3-momcntum cutoff in the NJL 
model. fj,Q is the quark number chemical potential. 

Given the thermodynamic potential, our job is to minimize it with respect to the fields tr, $ and $, and calculate 
the necessary quantities with the values of the fields so obtained. 



III. FORMALISM 



Here we give the details of the methods by which we have obtained the various thermodynamic quantities from the 
PNJL model. 



A. Taylor expansion of Pressure 



The pressure as a function of temperature T and baryon chemical potential fxo is given by. 



(4) 
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The first derivative of pressure with respect to /xq gives the quark number density. The second derivative is the 
quark number susceptibihty (QNS), which should show a power law divergence close to the CEP. On the Lattice 
where direct simulation for non-zero /io is not possible, the QNS and higher order derivatives (HODs) computed at 
fiQ = are used as Taylor expansion coefficients to extract chemical potential dependence of pressure. In fact the 
convergence of such an expansion has been tested with the HODs to obtain the CEP. At this point we should mention 
that in the present model the isospin chemical potential fii has not been included, which will be necessary to study 
the isospin asymmetric QCD matter. Also, the presence of diquark physics needs to be incorporated in the model to 
have a complete study of PNJL model for large /io- 

So, given as in the PNJL model we first solve numerically for a, <i> and $ using the following set of equations. 
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The values of the fields so obtained can then be used to evaluate all the thermodynamic quantities in mean-field 
approximation. To obtain the transition point we need to look at the behaviour of the temperature derivatives of the 
fields. The method used here to obtain the quantities d^/dT and da/dT is as follows. First we numerically solve 
Eqn. (jSJ for each value of temperature T for /ip — 0. The temperature difference between consecutive data is 0.1 
MeV. Then we obtain the slope of the (T, $(T)) and the {T,a{T)) curves. However, instead of using the difference 
method of extracting the slope, we do a fit to a Taylor expansion of the fields as a function of T, around the point 
where the slopes of the fields are required. We used a quadratic fitting function and the first order coefficient then 
gives us the required slope. We did the analysis more carefully near the transition where the slopes were obtained for 
every 1 MeV difference. Our results are identical to that obtained in Ref . |6^ . 

The field values obtain from Eqn.lO are then put back into to obtain pressure from We can then expand 
the scaled pressure as, 



rp4 



n=0 



(6) 



where. 



1 a" {PiT,fio)/T^ 



9 m 



Mo=0 



(7) 



We shall use the expansion around /io = 0. In this expansion, the odd terms vanish due to CP symmetry. We 
extract the expansion coefficients upto eight order. This has been motivated by the fact that on the Lattice the most 
recent results are also obtained upto this order. 

In general, to obtain the Taylor coefficients of pressure, one can use either of the two methods: 
Method (a): First the pressure is obtained as a function of fio for each value of T, and then fitted to a polynomial in 
liQ. The quark number susceptibility (QNS) and all other higher order derivatives (HODs) are then obtained from 
the coefficients of the polynomial extracted from the fit. 

Method (6): First obtaining the expressions for the derivatives of the pressure with respect to /io for the Taylor 

coefficients and then use the values of cr, and $ at zero chemical potential into these expressions. 

In any exact computation, these two methods should yield identical results. The Lattice however at present cannot 
use method (a) due to the complex determinant problem. On the other hand since we are using the mean field 
analysis, method (6) would give us wrong results as the mean fields used would be insensitive to /io. In this work we 
have computed all the observables using method (a). We have expanded the pressure in an eighth order polynomial 
in /io with even terms only (the odd terms should be zero due to CP symmetry). 

B. Cv, and C 

Given the thermodynamic potential fi, the energy density e is obtained from the relation. 



djn/T) 
dT 



= -T 



df 



(8) 
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The rate of change of energy density e with temperature at constant volume is the specific heat Cy which is given 
as, 



Cv = 



df 



= - T 



gf2 



(9) 



For a continuous phase transition one expects a divergence in Cy, which, as discussed earher, wiU translate into 
highly enhanced transverse momentum fluctuations or highly suppressed temperature fluctuations if the dynamics in 
relativistic heavy-ion collisions is such that the system passes close to the CEP. 

The square of velocity of sound at constant entropy S is given by. 



dP 



dP 
df 



df 



dn 

df 



T 



^dHl 

Qf2 



(10) 



Since the denominator is nothing but the C'v , a divergence in specific heat would mean the velocity of sound going 
to zero at the CEP. 

The conformal measure C is given by. 



C = — ^ ~ 1 - . 11 

e e 

Thus, a minima in the velocity of sound as expected near a phase transition or crossover may translate to a maxima 
of C. Considering the last relation we see that at asymptotic temperatures where the goes to the ideal gas value 
of 1/3, the conformal measure should go to the conformal limit C — 0. 

Given the relations H10|l and (|lll) we perform the same exercise of obtaining the n(T, /io = 0) from the PNJL 
model. We then obtain the Taylor expansion coefficient of 0(T, /^o = 0) around the temperatures at which we wish 
to obtain the values of Cy, v"^ and C. The expansion is done upto second order and fitted with the numerical values. 
From the first and second coefficients of this expansion we obtain dil/dT and d^fl/dT^. These are then put into the 
corresponding expressions for the thermodynamic quantities. 



IV. RESULTS 



As mentioned earlier we use the parameterization of the PNJL model as given in Ref.|63 and reproduced the 
behaviour of d^/dT and da/dT as shown in Fig^ For /ig = 0, the transition is quite sharp at Tchirai — 229 MeV 
and Tdeconfinement — 225 McV. For our subsequent results we shall present the quantities as a function of temperature 
in units of a cross-over temperature which is taken to be = 227 MeV. In view of the limitations of the model at 
high temperatures as mentioned earlier, we obtain results for thermodynamic quantities upto about 2.5Tc. 



A. Taylor expansion of Pressure 



We now present the pressure and its derivatives with respect to at /xq = 0. The pressure was fitted to a 
polynomial in /ip using the GNU plot program at different values of temperature. The maximum range of /xq was 
chosen to be 200 MeV. However near Tc the (which in this case is same as the least square) of the fit varies rapidly 
with the variation of range of over which the fit was done, and the actual range was chosen from minima of . 
The data points were spaced by 0.1 MeV. 

In table HI we present the fitted values of the Taylor coefRcients for a few values of temperature to show the 
dependence of these fitted coefficients on the range of and number of terms in the polynomial. Due to limitations 
of numerical accuracy and time, we chose to take a maximum of only upto the eighth order term in /iq. We hope to 
obtain higher order coefficients in future. The coefficients seem to be quite robust. For each temperature in the table 
the coefficients in the first case was actually used in the final result. 

Fig 13 displays the pressure scaled with that of Stefan-Boltzmann (SB) gas [P/Psb) as a function of T/Tc- Upto 
2.hTc the pressure grows from almost zero at low temperatures to about 90% of its ideal gas value. This is a bit high 
when compared to continuum estimates on the Lattice [tOI , which is about 80% of Psb ■ A comparison of our results 
with the PNJL result in Fig. 7a of Ref.|63|, shows a near perfect match, though Tq = 270 MeV in our measurements 
compared to their value of Tq = 190 MeV. 
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FIG. 1: d^/dT and da/dT as a function of temperature T (GeV) (also obtained in [6^ Fig.4b).The symbols show our data, 
the lines are spline fitting. 
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TABLE L Representative values of fitted Taylor coefficients to show the quality of fit. For each temperature, data in first line 
is used for the figures. (Here is same as least square.) 

Next we show the comparison of the coefficients C2, C4 and Cg obtained by us with Lattice data available in Table 
3.2 of Ref.j2Si|. FigOlshows the variation of the QNS C2 with T/T^- This shows an order parameter-like behaviour. 
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FIG. 2: Pressure scaled with Stefan-Boltzmann pressure as function of temperature for = from PNJL model. Arrow on 
the right indicates ideal gas value. 
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FIG. 3: The QNS as a function of T/Tc- Symbols are Lattice data [2^ . Arrow on the right indicates ideal gas value. 

Similar behaviour has been observed in another model study [t^ using the Density Dependent Quark Mass (DDQM) 
model. At higher temperatures the C2 reaches almost 85 % of its ideal gas value, consistent with Lattice data. 

The fourth order derivative C4, which can then be thought of as the "susceptibility" of C2 shows a peak aX T — T^ 
(Fig01). Near the transition temperature Tc, the effective model should work well and we observe that the structure 
of C4 is quite consistent with present day Lattice data H^ll^. Just above Tc however, there is a significant difference 
between our results of C4 and that of Ref.^^. While the Lattice values converge to the SB limit, ours is almost double 
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FIG. 4: The C4 as a function of T /Tc. The solid hne is from PNJL model, Symbols are Lattice data j2^. Arrow on the right 
indicates ideal gas value. 
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FIG. 5: The ce as function of T/Tc. The solid line is from PNJL model, Symbols are Lattice data p^ . 

of that value and shows only a weak convergence towards the SB limit. However, our results are consistent with the 
values of Xio in Fig. 3 in Ref.[2^, though it has been argued that these data would come down to the SB limit 
once a correct continuum limit is taken. Note that in the SB limit both C2 and C4 have only fermionic contributions. 
We expect that because the coupling strength is still large in this temperature regime it is unlikely that C4 should 
go to the SB limit within T < 2.5Tc. Moreover, the quark masses used in Ref.j2^ is considerably large (m/T — 0.4) 
to expect fermionic observables to go to the SB limit. However, it is possible that our overestimation is due to the 
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FIG. 6: The cs as function of T/Tc from PNJL model. 



use of mean field approximation. Recently, in a quasi-particle model |72j the Taylor coefficients have been found to 
match the Lattice data quite accurately. This model uses C2 to fix the T and /x dependence of an effective coupling 
in the self energy and then predicts the values of the higher order coefficients. In contrast, the PNJL model uses only 
pure gauge results on the lattice to fix the parameters of the Polyakov loop effective potential and then predicts all 
the coefficients including C2. 

The HODs Cg and Cg are plotted in FigOand Fig respectively. At high temperatures both of these HODs converge 
to zero. The structure near Tc is more interesting. Compared to other models the pattern of the HODs show much 
better resemblance to the Lattice results. The hadron resonance gas model describes the Lattice data well below Tc, 
but fails for T > Tc [tJ ■ The recently proposed |7^J scenario of coloured bound states also compares the HODs from 
Lattice data. However, their comparison at the present state is still not completely satisfactory. 



We now discuss the results for the other thermodynamic quantities that we have computed. First we consider the 
specific heat. As shown in FigQ Cy grows with increasing temperature and reaches a peak at Tc- Then it decreases 
sharply for a short range of temperature. Thereafter it shows a broad but shallow bump around T = 270 MeV and 
then gradually goes to a value little below the ideal gas value. Comparing with the recent lattice data for pure glue 
theory (Fig.Sf, Ref.js^), we see that the difference from the ideal gas value is far less in our case. For comparison, 
we have also plotted the values of 4e/T'*, at which the specific heat is expected to coincide for a conformal gas. As is 
clear from the plot, the values almost coincide at the large temperatures but not as perfectly as in Ref . [56| . 

Let us now consider the speed of sound and the conformal measure. In view of the possible relation between the 
two as indicated in we plot both these quantities together in Fig|Sl As can be seen from the figure, the two 

quantities indeed behave in such a correlated manner. We note here that C has been computed from the zeroth and 
first order coefficients of i7(r, /io = 0), whereas the has been measured from the first and second order coefficients. 
However the relation implied is not exact as can be seen in the figure. The value of P/e matches with that of for 
T < Tc and also goes close again above 2.5Tc. But in between these two limits the v'^ is distinctly greater than P/e. 
Thus, C would go to zero much faster if we replace P/e by v'^ for computing C. 

The v'^ is close to its ideal gas value at the temperature of about 2.5Tc. This is close to the results for pure glue 
theory on the Lattice as reported in Ref.l5a| and also that with 2 flavour Wilson Fermions in Ref. 75], with 2+1 
flavours of staggered quarks reported in Ref. |7^ , and with improved 2 flavour staggered fermions (P/ e was measured 
in this case). This possibly refers to an intriguing fact that is dominated by the gluonic degrees of freedom at least 



B. Cv, v1 and C 
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FIG. 7: Cv/T^ and 4e/r'' as function of T/Tc The arrow on the right shows the ideal gas value. 
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FIG. 8: Squared velocity of sound v1 and conformal measure C — A/e as function of T/Tc- The arrow on the right shows the 
ideal gas value for v^. For comparison with we also plot the ratio P/e. 



for large temperatures. However, near Tc the in Ref.|53 goes to a minimum value above 0.15, whereas we find the 
minima going close to 0.08, consistent with simulations with dynamical quark in Refs. pTSllT^ . and remarkably close 
to the softest point P/e = 0.075 in Ref. {IT\ . This is expected as the scalar a is supposed to be important in chiral 
dynamics for small quark masses. Even upto temperatures as low as 0.5rc the v1 in our case doesn't go much above 
0.1. This is also in contrast with the values reported in Ref. 54], where a confinement model [t^ has been used. These 
authors find a value of about v1 — 0.15 near and a value of 0.2 near 0.5rc- On the other hand v1 (isothermal) 
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measured in Ref . |7l| using the DDQM model, shows a better agreement with our results. 

The values for the conformal measure C also closely resembles the Lattice data of Ref. [s^l • Near T^, there is a slight 
departure. We find the value ranging from 0.7 to 0.75, whereas those authors find it to vary around 0.6 to 0.7 for 
pure glue theory. However, note that the dip at temperatures less than Tc is prominent in both the cases. At even 
lower temperatures we find C to increase. For a non-relativistic ideal gas, the ratio of P/e should go to zero, and thus 
C should then go to 1. Thus at lower temperatures such an increase is indeed expected. On the other hand at high 
temperatures either an ideal gas or a conformal behaviour should be recovered for which C should go to zero. Our 
results in this respect resemble the pure gauge Lattice results very closely. 

V. DISCUSSIONS AND SUMMARY 

We have studied several thermodynamic quantities of recent interest in QCD, in the framework of Polyakov-loop 
extended NJL (PNJL) model. Inclusion of Polyakov loop in the NJL model covers both the confinement as well as 
chiral symmetry breaking; the two most important characteristics of low energy QCD. The Polyakov-loop provides 
with the order parameter for deconfinement transition while the quark condensate acts as the order parameter for 
chiral transition. 

It is well known that the measurement of various fluctuations is of crucial importance to understand the characteris- 
tics of transition of confined matter to deconfined state. So, using the PNJL model, we have computed the derivatives 
of pressure upto eight order. It is expected that the susceptibilities might provide the most direct evidence for the 
order of QCD phase transition. In fact, a well defined peak in the susceptibility can indicate the cross over transition. 
On the other hand, a sharp diverging behaviour would indicate the existence of CEP. 

In the present paper, we have kept our computation close to /io = axis. The second derivative of pressure C2 rises 
steeply near Tc and then saturates near T = 2Tc. This behaviour indicates that near the transition the quark number 
fluctuations increase with temperature and then gradually saturates at higher temperatures. Such a behaviour is 
expected both from lattice as well as phenomenological studies. C4 has a peak around Tc, whereas, both cg and Cg 
show rapid variation across Tc- Considering the inherent differences of model and lattice approaches, our results 
are in good agreement with the lattice data. There is however a significant difference found in the values of C4 in 
the temperature range starting from a little above Tc- We speculate this may have occurred due to the neglect of 
fluctuations around the mean field in our computations; though one should keep in mind the very large masses used 
in Lattice computations in contrast to the value of 5.5MeV used in the PNJL model. 

We have also calculated the specific heat, speed of sound and the conformal measure. These results again are in 
close agreement with Lattice data. There seems to be a dominance of gluonic contribution at high temperatures. 
Near Tc however there are small but significant difference in Lattice data with pure QCD and that with dynamical 
quarks. Our results in this temperature region match well with the latter set. 

There are several issues that need to be addressed for a better understanding of the physics of QCD phase transition. 
NJL model works with constant four point coupling strength, the running coupling being averaged over a limited 
low energy kinematic domain. To make contact with QCD at high temperatures this coupling should have proper 
temperature dependence. Moreover, the parameters in the Polyakov-loop part should be tuned to get both the 
deconfinement and chiral transition within a temperature difference of couple of MeVs at most and also below, say, 
200 MeV, as indicated by the Lattice data. Note that similar parameter fitting of the Polyakov-loop in the PNJL 
model led the authors in Ref.js^ to couple the Polyakov-loop and SU{3) x SU{3) linear sigma model to study the 
K/tt ratio. However, there are more predictions in the present model, the most significant being the possibility of 
studying the QCD transition at finite baryon densities. Above and all one should also worry about the mean-field 
approximation involved here. For the Polyakov loop part, it would be worthwhile to explore approches similar to 
Ref.^3, where a local Polyakov loop has been coupled to the chiral Lagrangian. 

We have not yet studied the system with large chemical potentials. As already noted in Ref. [63l| . the PNJL model 
does have the requisite physics of a first order transition for low temperatures and high chemical potential, a cross-over 
for /io = and thus a CEP for some T = Te and /io — (J-Oe ■ The details in this case may be compared to other 
NJL model studies with high ^q. For example the chiral symmetry broken phase and normal quark matter phase of 
Ref. 79] show a similar value of the CEP. Moreover it has been shown in Ref. that with increase of strength of the 
diquark coupling the CEP moves towards higher temperatures and slightly lower values of fxo- Hence, the inclusion 
of diquark physics may be important to assess the exact location of the phase boundary as well as the CEP. 

A natural extension of our work would be to include the isospin chemical potential and study the behaviour for 
isospin and electric charge fluctuations. We hope to address some of these issues in future. 
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